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Abstract 

In  this  paper,  a general  technique  for  evaluation  of  measurements  by  the  method  of 
Least  Squares  is  presented.  The  input  to  the  method  consist  of  estimates  and  associated 
uncertainties  of  the  values  of  measured  quantities  together  with  specified  constraints 
between  the  measured  quantities  and  any  additional  quantities  for  which  no  information 
about  their  values  are  known  a priori.  The  output  of  the  method  consist  of  estimates  of 
both  groups  of  quantities  that  satisfy  the  imposed  constraints  and  the  uncertainties  of 
these  estimates.  Techniques  for  testing  the  consistency  between  the  estimates  obtained  by 
measurement  and  the  imposed  constraints  are  presented.  It  is  shown  that  linear  regression 
is  just  a special  case  of  the  method.  It  is  also  demonstrated  that  the  procedure  for 
evaluation  of  measurement  uncertainty  that  is  currently  agreed  within  the  metrology 
community  can  be  considered  as  another  special  case  in  which  no  redundant  information 
is  available.  The  practical  applicability  of  the  method  is  demonstrated  by  two  examples. 


1 Introduction 

In  1787,  the  French  mathematician  and  physicist  Laplace  (1749-1827)  used  the  method 
of  Least  Squares  to  estimate  8 unknown  orbital  parameters  from  75  discrepant  observa- 
tions of  the  position  of  Jupiter  and  Saturn  taken  over  the  period  1582-1745.  Since  then, 
the  method  of  Least  Squares  has  been  used  extensively  in  data  analysis.  Like  Laplace, 
most  people  use  a special  case  of  the  method,  known  as  unweighted  linear  regression.  The 
calculation  of  the  average  and  the  standard  deviation  of  a repeated  set  of  observations 
is  the  most  simple  example  of  that.  The  unweighted  regression  analysis  is  based  on  the 
assumptions  that  the  observations  are  independent  and  have  the  same  (unknown)  vari- 
ance. In  addition,  the  linear  regression  is  based  on  the  assumption  that  the  observations 
can  be  modelled  by  a function  that  is  linear  in  the  unknown  quantities  to  be  determined 
by  the  regression  analysis.  For  most  measurements  carried  out  in  practice,  none  of  these 
assumptions  can  be  justified.  In  order  to  evaluate  the  result  of  a general  measurement, 
in  which  some  redundant  information  has  been  obtained,  one  therefore  has  to  Apply  the 
method  of  Least  Squares  in  its  general  form. 

This  paper  describes  how  measurements  can  be  evaluated  by  the  method  of  Least 
Squares  in  general.  The  paper  is  based  on  an  earlier  work  of  the  author  [2]  but  includes 

1 Address:  Building  307,  Matematiktorvet,  DK-2800  Lyngby,  Denmark 
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several  new  features  not  published  before  as  well  as  practical  examples  from  the  daily 
work  at  DFM.  An  alternative  approach  is  described  in  [6]. 

2 Measurement  model 

In  a general  measurement,  a number  m > 0 of  quantities  is  either  measured  directly  using 
measuring  instruments  or  known  a priori,  for  example  from  tables  of  physical  constants 
etc.  The  (exact)  values  of  these  m quantities  are  denoted  Q 

C = 

Due  to  measurement  uncertainty,  the  values  z obtained  by  the  measurement  (or  from 
tables  etc.) 

Z — (-^l,  *•.,  2m) 

are  only  estimates  of  the  values  £.  The  standard  uncertainties  of  the  estimates  Zi, 

u(zi)  , i = 1, . . . ,m, 

are  determined  in  accordance  with  the  GUM  [1]  and  depend  on  the  accuracy  of  the 
instruments  and  the  reliability  of  any  tabulated  value  used.  In  general,  some  of  the 
estimates  Zi  may  be  correlated.  If  r(zi,Zj)  is  the  correlation  coefficient  between  the 
estimates  2;  and  zj  then  the  covariance  u(zi,zj)  between  these  two  estimates  is  given  by 


u{zi,zj)  =u(zi)r(zi,zj)u(zj). 

Because  of  the  uncertainty,  the  estimates  z can  be  considered  as  an  outcome  of  a Tri- 
dimensional random  variable  Z with  expectation  £ (the  exact  values  of  the  quantities) 
and  covariance  matrix  £ 


£ = u(  z,  z 


T\  _ 


/ u2(zx) 
u(z2,Zx) 


u(z1,z2) 

U2{z2) 


u(zi,zm)  \ 
u{z2,zm) 

U2(zm ) J 


In  addition  to  the  m quantities  for  which  prior  information  is  available  either  from 
direct  measurement  or  from  other  sources,  a general  measurement  may  involve  a number 
k > 0 of  quantities  for  which  no  prior  information  is  available.  The  values  of  these 
quantities  are  denoted  by 

/3  = (/?!,-  • - ,Pk)T  ■ 

In  general,  the  values  /3  and  £ are  constrained  by  a number  n of  physical  or  empirical 
laws.  These  constraints  may  be  written  in  terms  of  an  n-dimensional  function 


( /i(M  \ 

( 0 \ 

/2(/3,C) 

0 

— 

: 

l Suit 3,0  ) 

V 0 J 

k < n < m + k. 


(2.1) 


It  is  assumed  that  /j  : Q — » R , i = 1 . . . n,  are  differentiable  functions  (with  con- 
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tinuous  derivatives)  defined  in  a region  fi  C Rk+m  around  (13,  Q.  As  indicated  in  (2.1), 
the  number  n has  to  be  larger  than  or  equal  to  the  number  k of  quantities  for  which  no 
prior  information  is  available;  otherwise  some  of  the  values  (3  cannot  be  determined.  In 
addition,  the  number  n of  constraints  has  to  be  smaller  than  the  total  number  H m of 
quantities  involved;  otherwise  the  values  of  f3  and  C would  be  uniquely  determined  by 
the  constraints  and  no  measurements  would  be  needed. 

The  estimates  z,  the  covariance  matrix  S and  the  n-dimensional  function  f(/3,  Q are 
the  input  to  the  general  Least  Squares  method.  It  should  be  stressed  that  no  probability 
distribution  has  to  be  assigned  to  the  input  estimates  z.  On  the  contrary,  if  a probability 
distribution  has  been  assigned  to  an  estimate,  it  should  be  used  to  calculate  the  mean 
value  and  the  variance  of  the  estimate  which  should  then  serve  as  input  to  the  Least 
Squares  method. 

Like  any  other  covariance  matrix,  the  covariance  matrix  u( z,zT)  = S is  positive 
semi-definite.  Otherwise,  at  least  one  linear  combination  xrz  of  the  estimates  z would 
have  negative  variance  u(xTz,zTx)  = xT£x.  In  the  following  it  is  assumed  that  E is 
positive  definite  and  therefore  non-singular. 


3 Normal  equations 

Least  Squares  estimates  /3  anfi  C °f  the  values  /3  and  £ are  found  by  minimizing  the 
chi-square  function 

x2(C;z)  = (z  - C)Ts-1(z  - C) 

subject  to  the  constraints 

f(/3,C)  = 0. 

It  is  convenient  to  solve  this  minimization  problem  by  using  Lagrange  multipliers  [5]: 
If  a solution  (/3,  £)  to  the  minimization  problem  exists,  the  solution  satisfies  the  equation 

V(Mli)f(/3,C,A;z)  = o 

where 

W C A;  z)  = (z  - CfS"1  (z  - 0 + 2ATf  0,  C) 

for  a particular  set  of  Lagrange  multipliers  A = (Ai, . . . , A„)T.  By  taking  the  gradient  of 
the  function  $,  the  following  n + m + k equations  in  (/3,  £,  A)  evolve: 


Vpf($,C)T\  = 0, 
— C)  + Vcf(/3,CTA  - 0, 
f(/3,C)  = 0, 


(3.1) 


where 


V/jf’ 


(Oh.  ...  Oh  \ 

1 dPx  dpk  ' 


I Sh.  ...  §h  , 
\ 60,  60k  ) 


and  Vrf 


/ 9h 
1 eci 


I Qfn 

\ aci 


\ 

acm  * 


Mil  , 
acm  / 


The  equations  (3.1)  are  called  the  normal  equations  of  the  Least  Squares  problem. 
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4 Solving  the  normal  equations 


If  (fli,Ci,  A,)  is  an  approximate  solution  to  the  normal  equations,  a refined  solution 
(/3i+1,£|+i,  Aj+i)  can  be  found  by  the  iteration 


A 

Ci 

o 


A/3,  \ 
A£, 
A,+ 1 / 


, Z = 1, ...  ,oo. 


The  step  (A/3,,  AC„  A,+i)  is  given  by 


/ A/3,  \ / 0(fcd)  \ 

D(/3„C,)  AC,  = IT^C  i)  , (4-1) 

V A,+1  ) \ -f(/3„C ,)  ) 

where 

/ 0(fc,fc)  V^f(/3„C,)T  \ 

D(/3„C,)=  0^)  5T1  Vcf(/3„C,)T  (4.2) 

V V^f(/3„C,)  Vcf(/3„C,)  ) 

is  a symmetric  matrix.  This  iteration  procedure  is  similar  to  Newton  iteration  except 
that  the  second  order  partial  derivatives  of  the  functions  /,  have  been  neglected  as  it  is 
practice  to  do  in  non-linear  Least  Squares  estimation  [4] . 

In  order  to  reduce  the  effects  of  numerical  rounding  errors,  it  is  recommended  to  calcu- 
late the  step  (A/3,,  AC,,  A,+i)  by  solving  the  linear  equations  (4.1)  by  Gauss-Jordan  elim- 
ination with  full  pivoting  [4],  This  algorithm  also  provides  the  inverse  matrix  D(/3,,  C,)_1 
which  is  needed  at  the  final  stage  for  estimating  the  covariance  matrix  of  the  solution  as 
shown  in  Section  5. 

If  proper  starting  values  (/31,Ci)  are  selected,  the  iteration  is  expected  to  converge 
towards  the  solution  (/3,  C) 


Since  the  solutions  C are  expected  to  be  close  to  the  estimates  z of  C available  a 
priori,  the  estimates  z are  obviously  the  proper  starting  values  Ci  to  be  selected  for 
the  iteration.  The  selection  of  proper  starting  values  (3l  is  more  difficult  in  general.  If, 
however,  f(/3,  C)  are  linear  functions  in  the  variables  /3,  the  iteration  process  will  converge 
after  a few  iterations,  independent  of  the  choice  of  /31 . 

Most  differentiable  functions  f(/3,  C)  can  be  handled  by  the  described  method.  In 
order  to  get  reliable  standard  uncertainties,  it  is  required,  however,  that  the  function 
can  be  approximated  by  a first  order  Taylor  expansion,  i.e. 

f(/3,  C)  = m C)  + V^f(/3,  C)(/3  -$)  + Vcf(/3,  C)(C  - C) 

when  the  values  /3  and  C are  varied  around  the  solution  /3  and  C on  a scale  comparable  to 
the  standard  uncertainties  of  the  solution.  If  this  vaguely  formulated  criterion  is  met,  the 
function  f(/3,  £)  is  said  to  be  linearizable.  Note  that  almost  any  differentiable  function 
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will  be  linearizable  if  the  standard  uncertainties  are  sufficiently  small.  On  the  other 
hand,  if  the  uncertainties  are  sufficiently  high,  all  non-linear  functions  will  no  longer  be 
linearizable.  The  requirement  that  f(/3,  C)  is  linearizable  is  considered  to  be  the  only 
major  limitation  of  the  method  of  Least  Squares! 

It  should  be  mentioned  that,  the  minimization  using  Lagrange  multipliers  will  fail  in 
case  the  gradients  Vpfc  and  of  one  of  the  constraint  functions  /,  are  both  equal  to 
zero  at  the  point  of  the  solution  (/3,  £),  This  gives  some  restrictions  on  how  a constraint 
may  be  formulated.  A function  /,  defining  a constraint  may,  for  example,  be  replaced  by 
the  square  of  that  function,  ff.  But  since  fi(j3,  Q = 0,  the  gradient  of  ff  will  be  zero  at 
the  point  of  the  solution  (/3,  £)  although  the  gradient  of  /,  is  not. 

5 Properties  of  the  solution 

Since  the  solution  (/3,£,  A)  depends  on  the  estimates  z,  which  are  considered  as  the 
realization  of  the  multivariate  random  variable  Z,  the  solution  (/3,£,  A)  can  also  be 
regarded  as  a multivariate  random  variable.  If  the  functions  /,  (/3,  Q are  linearizable,  the 
estimates  (/ 3 , C)  are  linear  in  Z 

( $ \ ( 0 \ ( 0<M)  \ 

l a / V o ) +D(AC)_1  ( =-£l70  j ' (51) 

In  that  case,  the  expectation  of  the  solution  is 


which  means  that  (/3,C)  are  central  estimators  of  the  values  (/3,C).  Under  the  same 
assumption,  the  covariances  of  the  solution  are  given  by  the  symmetric  matrix  D(/3,  C)-1 
provided  by  the  Gauss-Jordan  elimination  algorithm2 

( m(/3,/3T)  u0,£T)  ()(fc,n)  \ 

u(CJ f)  u(  C,CT)  0(m’n)  =D(/3,C)-1-D(AC)-1.  (5.2) 

V 0(n,fc)  ()(n,m)  -u(A,Ar) 

This  relation  can  be  derived  as  follows.  Partition  the  symmetric  matrix  D“'  into  nine 
sub-matrices  similar  to  the  left  hand  side  of  (5.2)  or  similar  to  the  partitioning  of  D 
according  to  the  definition  (4.2).  Express  the  covariance  matrix  of  the  solution  (5.1)  in 
terms  of  the  covariance  matrix  S of  the  random  variable  Z and  the  matrix  D-1 . Insert 
the  partitioned  D_1  into  the  resulting  matrix  double  product  and  express  the  covariances 
of  the  solution  in  terms  of  £~x  and  the  sub-matrices  of  D1.  Reduce  these  expressions 
to  the  final  result  by  multiple  use  of  the  nine  relations  between  the  sub-matrices  of  D1 
and  D derived  from  the  identity  DD1  = I. 

2The  empty  brackets  in  the  left  hand  matrix  indicates  the  parts  of  D-1  that  do  not  contain  information  about 
covariances. 
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From  equation  (5.1)  and  (5.2)  the  covariance  matrices  between  (0,  £)  and  the  estim- 
ates z are  found  to  be 

u(z,(3T)=u(t,$T) 
u(z,<T)  =u(C,CT). 

From  the  last  of  these  two  relations,  a relation  of  particular  interest  is  derived, 

u(z  - C,  zT  - f)  = zT)  - u(C,  <f). 

For  the  diagonal  elements,  this  relation  reads 

u2{zi-Ci)  = u2(zi)-u2(Q)  , i = 1, . . . ,m. 

That  is,  the  variance  of  the  difference  between  the  initial  estimate  Zi  of  Q and  the  refined 
estimate  Q is  equal  to  the  variance  of  minus  the  variance  of  Q.  This  relation  is  useful 
when  testing  if  the  difference  — Q is  significantly  different  from  its  zero  expectation. 

6 x2  test  for  consistency 

When  the  estimates  (/3,C)  have  been  found,  the  minimum  y2  value 

x2(C;z)  = (z-C)rs_1(z-c) 

can  be  used  to  test  if  the  measured  values  z are  consistent  with  the  measurement  model 
(2.1)  within  the  uncertainties  defined  by  the  covariance  matrix  E.  If  the  model  is  lin- 
earizable,  the  expectation  of  the  random  variable  x2(C;  Z)  is  equal  to  the  number  m of 
measured  quantities,  minus  the  number  m + k of  adjusted  quantities,  plus  the  number 
n of  constraints,  that  is 

E x2(C; Z)  =m—(m  + k)  + n = n-k  — u. 

If,  in  addition,  the  random  variables  Z are  assumed  to  follow  a multivariate  normal 
distribution  with  mean  values  £ and  covariance  matrix  E,  the  random  variable  x2(Ci  Z) 
will  follow  a x2M  distribution  with  u = n — k degrees  of  freedom.  In  that  case,  the 
probability  p of  finding  a x2  value  larger  than  the  value  x2(Ciz)  actually  observed  can 
be  calculated  from  the  x2(^)  distribution 

P = P[x2{v)  >X2(C,z)}  = 1 ~ P{x2{v)  <X2(C,Z)}- 

If  this  probability  p is  smaller  than  a certain  value  a,  the  hypothesis  that  the  meas- 
ured values  are  consistent  with  the  measurement  model  has  to  be  rejected  at  a level  of 
significance  equal  to  a.  As  the  result  of  measurements  are  normally  quoted  at  a 95%  level 
of  confidence,  an  a = 5%  level  of  significance  is  a reasonable  choice  for  the  consistency 
test. 

Although  the  assumption  of  a normal  distribution  of  Z may  not  be  fulfilled,  it  is 
suggested  to  carry  out  the  test  of  consistency  as  described  above  anyway.  This  is  justified 
by  the  fact  that  a value  of  X2(C;Z)  significantly  higher  than  the  expectation  v indicates 
inconsistency  no  matter  what  the  distribution  of  Z might  be.  The  calculated  probability 
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p simply  describes  how  unlikely  the  observed  y2  value  is  if  a normal  distribution  is 
assigned  to  Z. 

7 Normalized  deviations 

If  the  test  described  in  the  previous  section  leads  to  a rejection  of  the  measurements, 
a tool  for  identifying  the  outlying  measurements  is  desirable.  A measured  value  Z;  is 
defined  as  an  outlier  if  the  difference  Zi  - Q is  significantly  different  from  zero  taking 
into  account  the  standard  uncertainty  u(z{  — (i)  of  that  difference.  This  leads  to  the 
introduction  of  the  normalized  deviation  d,  defined  by3 

d,=  = Zi~Ci  = , i = 1, . . . ,m. 

The  normalized  deviation  d{  has  zero  expectation  and  variance  1.  A normalized  de- 
viation with  \di\  larger  than  2 or  3 is  therefore  rather  unlikely  no  matter  what  the 
distribution  of  the  random  variable  di  might  be. 

If  a multivariate  normal  distribution  is  assigned  to  Z and  the  model  function  f(/3,  Q 
is  linearizable,  the  normalized  deviation  di  is  normally  distributed, 

di  S N( 0, 1)  , i = 1, . . . ,m. 

In  that  case 

P {\di\  > 2}  = 5%, 

and  a measurement  with  |dj|  > 2 is  therefore  identified  as  an  outlier  at  a.  5%  level  of 
significance.  It  is  suggested  to  use  the  criteria  |d*|  >2  to  identify  potential  outliers  even 
if  the  distribution  assigned  to  Z is  not  normal. 

8 Adjustment  of  a variance  <r2 

If  some  values  Zi  have  a common  but  unknown  variance  u2{zf)  = cr2,  this  variance  can 
be  estimated  by  adjusting  cr2  by  an  iterative  procedure  until  the  “observed”  y2  value 
becomes  equal  to  its  expectation  value  v 

x2(C;  z)  = (z  - C)TS_1(z  - C)  = v, 

where  the  covariance  matrix  £ is  a function  of  the  unknown  variance  a2.  As  the  estimates 
C depends  on  the  value  assigned  to  cr2,  these  estimates  have  to  be  updated  together  with 
the  estimates  £ each  time  the  value  of  cr2  is  changed  during  the  iteration. 

This  way  of  estimating  the  unknown  variance  a2  leads  to  the  well-known  expression 
for  the  standard  deviation  in  the  case  of  a repeated  measurement  of  a single  quantity  as 
shown  in  Section  13. 


3If  u(zi  — fy)  = 0,  the  difference  Zi  — Q will  be  zero  as  well  and  d;  may  be  set  equal  to  zero.  This  situation  occurs 
whenever  there  is  no 'redundant  information  available  regarding  the  value  of  the  quantity  Cj. 
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9 Example:  Calibration  of  an  analytical  balance 

An  analytical  balance  with  capacity  Max= 220  g,  resolution  <2=0.1  mg,  and  built-in  ad- 
justment weight  was  calibrated  by  DFM  in  October  1999  during  an  inter-laboratory 
measurement  comparison  piloted  by  DFM.  Two  mass  standards  were  used  as  reference 
standards.  One  of  them  was  a traditional  200  g weight  (named  R200g)  of  known  con- 
ventional mass  value4  mu  and  density  pR.  The  other  reference  standard  was  a specially 
designed  200  g stack  of  weights  consisting  of  four  discs  (named  100g,  50g,  25g  and  25g*) 
machined  from  the  same  metal  bar  of  known  density  p.  The  conventional  mass  val- 
ues mi,  m2,  m3,  m 4 respectively  of  these  four  discs  were  not  known  a priori;  only  the 
conventional  mass  value  ms  = mi  + m2  + m3  + of  the  stack  was  known. 

The  calibration  was  performed  by  placing  a weight  combination  at  the  weighing  pan 
of  the  balance  and  by  recording  the  corresponding  average  indication  I in  the  display. 
A total  of  18  weight  combinations  were  used.  Each  weight  combination  was  weighed  3 
times  from  which  the  average  indication  was  calculated.  The  calibration  was  repeated  4 
times  during  a period  of  10  days  in  which  the  inter-laboratory  comparison  took  place. 
From  these  four  calibrations,  a grand  average  indication  Ii,i  = 1, . . . , 18  was  calculated 
for  each  of  the  18  weight  combinations  specified  in  Table  1.  The  standard  uncertainty  of 
the  grand  average  was  estimated  from  the  observed  variation  in  indication  over  the  four 
calibrations. 
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25g 

25g 

25g 

25g* 

25g* 

25g* 

25g* 

ho 

In 

1 12 

I\Z 

/l4 

hb 

Il6 

Il7 

ha 

50g 

50g 

50g 

50g 

25g 

25g 

25g* 

R200g 

R200g 

25g 

25g 

25g* 

25g* 

25g* 

Tab.  1.  The  weight  combinations  corresponding  to  the  18  balance  indications 


Due  to  the  effect  of  air  buoyancy,  the  balance  indication  depends  not  only  on  the 
mass  of  the  weighed  body,  but  also  on  the  density  of  the  body  as  well  as  the  density  of 
the  air.  When  calibrated  in  air  of  known  density  a,  the  reference  indication  Ir  of  the 
balance  corresponding  to  a load  generated  by  a weight  with  conventional  mass  value  m 
and  density  p is  given  by 

Ir  = m ( 1 - (a  - a0)  ( - - — 

\ \P  Po 

where  oo=1.2  kg/m3  and  po=8000  kg/m3  are  the  reference  densities  of  the  air  and  the 
weight  respectively  to  which  the  conventional  value  of  mass  refers.  As  a model  for  the 


4The  conventional  mass  value  of  a body  is  defined  as  the  mass  of  a hypothetical  weight  of  density  8000  kg/ 
that  balances  the  body  when  weighed  in  air  of  density  1.2  kg/m3  and  temperature  20  °C. 
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ms 

[g] 

mn 

[g] 

PR 

[kg/m3] 

[kg/m3] 

[kg/m3] 

h 

[div] 

z 

199.988816 

199.999043 

7833.01 

7965.76 

1.1950 

199.988617 

u{z) 

0.000010 

0.000008 

0.29 

0.71 

0.0035 

0.000023 

c 

199.988814 

199.999043 

7833.01 

7965.76 

1.1946 

199.988620 

w(C) 

0.000010 

0.000008 

0.29 

0.71 

0.0035 

0.000011 

d 

1.66 

-1.66 

1.66 

-1.66 

1.66 

-0.16 

h 

h 

h 

h 

h 

It 

[div] 

[div] 

[div] 

[div] 

[div] 

[div] 

z 

199.988608 

174.992133 

175.009992 

150.013558 

149.980675 

125.002083 

u{z) 

0.000023 

0.000023 

0.000023 

0.000023 

0.000023 

0.000023 

c 

199.988620 

174.992149 

175.010024 

150.013558 

149.980672 

125.002087 

u(0 

0.000011 

0.000012 

0.000012 

0.000013 

0.000012 

0.000014 

d 

-0.56 

-0.77 

-1.61 

0.03 

0.14 

-0.20 

h 

h 

ho 

In 

1\2 

il3 

[div] 

[div] 

[div] 

[div] 

[div] 

[div] 

z 

124.984217 

100.005650 

99.982925 

74.986433 

75.004325 

50.007892 

u(z) 

0.000023 

0.000023 

0.000023 

0.000023 

0.000023 

0.000023 

c 

124.984212 

100.005632 

99.982899 

74.986450 

75.004325 

50.007881 

«(<) 

0.000014 

0.000013 

0.000013 

0.000014 

0.000014 

0.000012 

d 

0.25 

0.93 

1.38 

-0.87 

0.03 

0.54 

/l4 

hs 

he 

Il7 

hs 

[div] 

[div] 

[div] 

[div] 

[div] 

z 

49.974992 

24.978533 

24.996417 

199.998867 

199.998875 

u{z) 

0.000023 

0.000023 

0.000023 

0.000023 

0.000023 

c 

49.974995 

, 24.978557 

24.996432 

199.998851 

199.998851 

«(<) 

0.000013 

0.000011 

0.000011 

0.000011 

0.000011 

d 

-0.19 

-1.17 

-0.77 

0.78 

1.19 

/ 

A 

m i 

7712 

m3 

7774 

[g/div] 

[1/div] 

tel 

fg] 

[g] 

fel 

0 

1.00000186 

-4.4E-09 

100.005774 

50.007963 

24.978601 

24.996476 

u(P) 

0.00000019 

1.0E-09 

0.000011 

0.000010 

0.000010 

0.000010 

Tab.  2.  Measured  and  estimated  values  and  associated  standard  uncertainties. 


calibration  curve  of  the  balance,  a second  order  polynomial  through  zero  is  assumed 

IR  = f{l  + AI 2) 

where  / and  A are  unknown  quantities  to  be  determined  from  the  calibration  data. 

In  this  example,  there  are  m — 23  quantities  for  which  prior  information  is  available 
from  the  measurements  performed: 

C = (ms,  fUR,  PRi  Pi  a,  I\ , • • • , hs)T 

whereas  there  are  k = 6 quantities  for  which  no  prior  information  is  available: 

/3  = (f,A,m1,m2,Tn3,mA)T. 
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/ 

A 

mi 

m2 

m3 

7714 

/ 

1 

-0.945 

0.021 

0.071 

0.096 

0.096 

A 

-0.945 

1 

0.124 

-0.016 

-0.094 

-0.094 

mi 

0.021 

0.124 

1 

-0.194 

-0.269 

-0.268 

m2 

0.071 

-0.016 

-0.194 

1 

-0.287 

-0.287 

m3 

0.096 

-0.094 

-0.269 

-0.287 

1 

-0.287 

m4 

0.096 

-0.094 

-0.268 

-0.287 

-0.287 

1 

Tab.  3.  Correlation  coefficients  of  the  estimated  $ values. 


Between  these  quantities,  there  are  n = 19  constraints: 

/ (mi  +m2  + m3  + m4)  (l  - (a  - «o)  -f(h  + All)  ^ 


f(/3,C) 


mR  (l  - (o  - no)  (hs  + AI2m) 

ms  - (mi  + ra3  + m3  + m4) 


= 0 . 


The  measured  values  z and  associated  standard  uncertainties  are  given  in  Table  2 
under  the  row  headings  2 and  u(z).  All  measured  values  are  assumed  to  be  uncorrelated. 

By  solving  the  normal  equations,  the  estimates  C and  /3  and  associated  standard 
uncertainties  given  in  Table  2 under  the  row  headings  (,  u((),  /3  and  u0)  are  obtained. 
Selected  correlation  coefficients  derived  from  D(/3,  C)-1  are  given  in  Table  3.  The  ob- 
served minimum  x2  value  is  x(C  z)  = 8.6  which  should  be  compared  to  the  expectation 
value  v = n — k = 19  — 6 = 13.  Since  P {x2(13)  > 8.6}  = 80.3%,  it  is  concluded  that 
the  measured  values  are  consistent  with  the  specified  constraints  taking  into  account  the 
measurement  uncertainties.  This  conclusion  is  confirmed  by  the  calculated  normalized 
deviations  given  in  Table  2 under  the  row  heading  d;  all  normalized  deviations  satisfy 
the  criterion  |d|  < 2. 

From  the  estimates  of  the  quantities  / and  A and  the  associated  covariance  matrix, 
the  error  of  indication  E,  defined  as 

E = I-IR  = I-f(l  + AI 2) , 


and  the  associated  standard  uncertainty  u(E)  can  be  calculated  as  a function  of  the 
indication  I.  The  result  is  shown  in  Figure  1 as  the  full  lines  representing  E — u(E ),  E, 
and  E+u(E).  The  measured  points  25,-,  i = 1, . . . , 18  shown  in  the  figure  are  the  observed 
average  balance  indications  J,  minus  the  corresponding  reference  values  Jr.  The  error 
bars  of  the  measured  points  indicate  the  standard  uncertainties  u(Ei)  that  have  been 
calculated  taking  into  account  the  covariance  between  I,  and  IR. 


10  Example:  Evaluation  of  calibration  history 

A weight  (named  Rlmg)  of  nominal  mass  1 mg  has  been  calibrated  39  times  in  the 
period  1992-2001.  For  calibration  number  i,  the  mass  m,  of  the  weight  at  the  time 
U and  the  associated  standard  uncertainties  u(rrii)  and  u(ti)  are  given.  The  calibration 
history  of  the  weight  is  shown  in  Figure  2 as  dots  with  error  bars  indicating  the  standard 
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Fig.  1.  Error  of  indication  of  the  calibrated  balance. 


uncertainties;  the  scale  mark  1992-01  on  the  time  axis  indicates  the  position  of  the  date 
1 January  1992  etc. 

Due  to  wear  and  changes  in  the  amount  of  dirt  adsorbed  to  the  surface,  the  mass  of 
the  weight  is  expected  to  change  in  time.  A reasonable  model  of  the  change  in  mass  as  a 
function  of  time  is  a superposition  of  a deterministic  linear  drift  and  a random  variation 

mi  = di  + atfi  + 6mt  , t = l,...,39, 

where  Jm*  is  a random  variable  with  zero  expectation  and  variance  cr2.  The  drift  para- 
meters di , <22  and  the  associated  covariance  matrix  as  well  as  the  variance  a2  are  unknown 
a priori  and  are  to  be  estimated  from  the  calibration  history  available.  Once  the  estim- 
ates di  and  &2  have  been  found,  it  is  possible  to  predict  a value  m of  the  mass  of  the 
weight  as  a function  of  time  t 

m = a i + a2t  + 6m , 

where  6m  = 0 with  standard  uncertainty  u(6m ) = a.  The  standard  uncertainty  of  the 
predicted  mass  value  is  given  by 

u2(rh ) = it2(di)  + t2u2(a2)  + 2tu{a\,  02)  + o2 . 

The  measurement  model  used  for  evaluating  the  calibration  history  is 

c = (mi,...,m39,ti,...,t39,5mi,...,<5m39)r  , (3  = {aua2)T , 


1992-01  1993-01  1994-01  1995-01  1996-01  1997-01  1998-01  1999-01  2000-01  2001-01 


Fig.  2.  Evaluation  of  the  calibration  history  of  a 1 mg  weight  assuming  that  a = 0. 


Fig.  3.  Evaluation  of  the  calibration  history  of  the  1 mg  weight  with  a adjusted 
to  0.092  /ig. 


mo  = 


mi  - (ai  + 02*1  + Smi) 


m.39  — (Oi  + 02*39  + <5^39) 


The  measured  values  z are  given  by  the  calibration  history,  except  for  the  values  of 
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8m,i,i  = 1, . . . , 39  which  are  set  equal  to  the  expectation  value  zero.  The  associated  cov- 
ariance matrix  u( z,zT)  = £ is  built  up  from  the  uncertainties  i/(m, ) and  u(t,)  available 
from  the  calibration  history  and  a negligible  but  finite5  initial  value  of  the  unknown 
variance  a2.  Since  the  standard  uncertainties  u(m,)  are  of  the  order  0.1  pg.  the  value 
<7=lE-07  pg  is  considered  negligible  and  is  selected  as  a starting  point. 

By  solving  the  normal  equations,  estimates  ai  and  o2  of  the  drift  parameters  and  the 
associated  covariance  matrix  are  found  after  a few  iterations.  The  predicted  value  m of 
the  mass  of  the  weight  and  the  associated  standard  uncertainty  u(m)  as  a function  of 
time  are  shown  in  Figure  2 as  solid  lines.  The  normalized  deviations  d associated  with 
the  mass  values  m*  are  shown  in  Figure  2 as  well6.  The  observed  minimum  chi-square 
value  is  x2  = 66-6  which  is  large  compared  to  the  expectation  value  v = 39  - 2 = 37. 

Since  P{x2(37)  > 66.6}  = 0.2%,  the  hypothesis  a = 0,  or  no  random  variation  in  the 
mass,  is  rejected  at  a 0.2%  level  of  significance. 

The  value  of  a is  therefore  increased  as  described  in  Section  8 until  the  calculated 
minimum  x2  value  becomes  equal  to  its  expectation  value  v = 37.  In  this  way  the 
standard  uncertainty  reflecting  the  random  variation  of  the  mass  of  the  weight  is  found 
to  be  cr=0.092  pg.  The  result  of  the  evaluation  of  the  calibration  history  after  adjustment 
of  c7  is  shown  in  Figure  3.  Note  the  significant  increase  in  the  standard  uncertainty  of 
the  predicted  value  of  the  mass  of  the  weight  and  the  decrease  in  the  absolute  value  of 
the  normalized  deviations  d. 

The  calibration  history  can  also  be  evaluated  by  an  iterative  technique  based  on  linear 
regression  [3].  The  results  obtained  are  identical  to  the  results  presented  in  this  section. 

11  Case  I:  Univariate  output  quantity,  Y = h(X i, . . . , X at) 

In  this  section  it  is  shown  that  the  evaluation  of  measurements  by  the  method  of  least 
squares  is  consistent  with  the  generally  accepted  principles  for  evaluating  measurement 
uncertainty  as  described  in  the  GUM  [1], 

Using  the  nomenclature  of  the  GUM,  a univariate  output  quantity  Y is  assumed  to 
be  related  to  N input  quantities  X\, . . . , through  a specified  function  h, 

Y = h(X1,...,XN). 

The  values  assigned  to  the  input  and  output  quantities  are  denoted  xi,...,xn  and 
by  y respectively.  In  the  nomenclature  of  this  paper,  the  measurement  model  is 

C = {Xi,...,Xn)t  , f3  = (Y), 
f(/3,C)  = (r-fi(X1....,AV))  = 0. 

The  measured  values  are 

z = (xu...,xn)t 

5If  the  variance  a2  is  assumed  to  be  exactly  zero,  the  quantities  <5m;  have  to  be  removed  from  the  model.  Otherwise 
the  covariance  matrix  £ will  be  singular. 

6The  absolute  value  of  normalized  deviations  of  t%  and  <5m.j  is  equal  to  the  absolute  value  of  the  normalized 
deviation  of  m;. 
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with  the  known  covariance  matrix 


/ u2(x i) 


S = u(z,zT)  — 

\ u(xN,xi) 

The  coefficient  matrix  D of  the  normal  equations  is 

0 O*1'*) 

D(/3,<)=(  O^’1)  IT1 


— Vxh(x)T 

1 -Vx/i(x)  0 

where 

„ , f dh  dh  \ 
x _ ■ 

In  the  present  case,  the  solution  to  the  normal  equations  is  found  after  one  iteration, 
y = $ = h(x  i,...,xN)  , C = {xi,...,xn)t,  A = 0. 


The  associated  covariances  are  given  by 


2(y)  u(y,C  ) 0(1,1)  \ 


u(Cy)  «( C,0  0W1) 

q(U)  q(1  ,N)  _u2(A)  J 


= D(/3,0 


-1 


Vxh(x)  X Vxh(x)T  Vxh(x)S  1 
SVxh(x)T  X 0 

1 0 0 


In  other  words, 
,2 


N N 


U 


(y)  = vxh(x)  II  Vx/i(x)T  = ^ Ciu(a:i, Xj)cj,  c,  = (x<) 

i=l  i=l 

which  is  identical  to  the  linear  variance  propagation  formula  given  in  the  GUM. 


12  Case  II:  Linear  regression,  Y = Xa 

Linear  regression  is  applied  when  there  is  a linear  relationship  Y = Xa  between  some 
observed  quantities  Y and  some  unknown  quantities  a.  The  design  matrix  X is  made  up 
of  known  elements  that  may  be  given  as  specified  functions  of  one  or  several  independent 
variables.  In  the  notation  of  this  paper,  the  measurement  model  for  the  linear  regression 
problem  is 

C = Y = (Yi,...,Y„)t  , /3  = a = (a1,...,afc)T, 
f«,/3)  = Y — Xa  = 0, 

where  X-n,fc)  is  the  known  design  matrix.  The  measured  values  are 

z = y = {yi,---,yn)T 
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with  known  covariance  matrix 

l u2{yi)  ■■ 

" u(yuyn)  \ 

E = u(  z,zT)  = : 

■ 

\ u(yn,y i)  ■■ 

1 • u2{yn)  j 

The  coefficient  matrix  D of  the  normal  equations  is 

/ Q{k,k)  Q(fc,n)  _XT  \ 

D(/3,<)  = I E-1  It”’”)  I . 


\ —X  l(">n)  ot"’”)  ) 

Again,  the  solution  to  the  normal  equations  is  found  after  one  iteration, 

a = /3  = CXrE-1y  , Y - C - Xa  , A = -E"1(Y-y) 
where  C = (XTE-1X)_1.  The  associated  covariances  are  given  by 
( u( a,ar)  u(a,YT)  Q(k,n)  \ 


u{  Y,Y?) 

()(n,n)  j 

= D(/9,C  )_1 

Q(n,n) 

-u(X,XT)  ) 
( C 

CXT 

-CXtE-1 

= 

xc 

XCXT 

I - XCXrE 

V -E_1XC 

I - e_1xcxt 

S_1XCXTE_1- 

that  is, 

a = CXrE-1y  , u(a,aT)  = C = (XTE-1X)-1 
as  is  known  from  the  theory  of  linear  regression. 

13  Case  III:  Repeated  observations  of  a single  quantity 

Assume  that  a quantity  X is  measured  n times  with  the  same  uncertainty  a.  Such  a 
measurement  can  be  modelled  by  n quantities  X\ , . . . , X„  having  a common  value  f i 

£ = X=:  (X1,...,Xn)T  , £ = (,*), 

/ Xi  - fi  \ 

mo=  =o. 

y Xn-  y j 

The  measured  values  are 

z = x = (xi,...,xn)T, 

and  under  the  assumption  that  the  measurement  results  are  mutually  independent,  the 
associated  covariance  matrix  is  given  by 

/ (T2  •••  0 \ 
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The  coefficient  matrix  D of  the  normal  equations  is 

(of1*1)  — i(h*0 

0(n,l)  ff-2j(n,n)  j(n,n) 

— iffi.1)  Q(n>") 

where  1 denotes  a matrix  with  all  elements  equal  to  1.  The  solution  of  the  normal 
equations  is  found  after  one  iteration, 

1 71 

fi  = 0=±Yxi  , X = C = £l(n,1)  , A = -ct-2(X  - x). 
n 1 ' 


i= 1 


The  associated  covariances  are  given  by 

1 «2(£)  w(£,XT)  ()(1,n)  N 

u(X,A)  u(X,  XT)  0(n’n) 

V 0(n,1)  0(n,n)  -u(\\T)  ) 


a2n  * 


D (AC)"1 


CT2n-li(i,») 


n-il(l,n) 


a2n-ll(n,l)  CT2n-li(n,n) 


j(n,n) 


■ n 


n 


llffid)  l(n.n)  _ n-ll(n.n)  0--2(?iJ-1i(n,n)  _ j(n,n)j 


As  expected, 


A=-Va:i  , u2{(i)=a2/n. 

n 

t=i 

If  cr2  is  not  known  a priori,  it  can  be  estimated  by  solving  the  equation 

*’<<;*)  = £k^  = n- 1, 


i= 1 


crz 


i.e., 


ct2  = ^31 S (**  - v)2 

i=  1 


which  is  the  well  known  expression  for  the  experimental  standard  deviation  s. 


14  Conclusion 

A general  technique  for  evaluation  of  measurements  by  the  method  of  Least  Squares 
has  been  presented.  The  applicability  of  the  method  has  been  demonstrated  by  two 
examples.  It  has  been  shown  that  the  method  is  fully  compatible  with  the  generally 
accepted  principles  for  evaluation  of  measurement  uncertainty  laid  down  in  the  GUM 
and  that  ordinary  linear  regression  is  just  a special  case  of  the  method. 

The  input  to  the  method  consists  of 

• An  estimate  of  the  value  of  each  measured  quantity,  including  any  relevant  influence 
quantity. 

• The  covariance  matrix  of  these  estimates  formed  by  the  standard  uncertainties  of 
the  estimates  and  the  correlation  coefficients  between  the  estimates. 
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• A measurement  model  describing  all  the  known  relations  between  the  measured 
quantities  and  some  additional  quantities  (if  needed)  for  which  no  prior  information 
is  available. 

The  output  of  the  method  consists  of 

• An  adjusted  estimate  of  the  value  of  each  measured  quantity  and  an  estimate  of 
each  additional  quantity  introduced  in  the  measurement  model. 

• The  covariance  matrix  of  all  these  estimates  from  which  the  standard  uncertainties 
and  correlation  coefficients  can  be  calculated. 

• A chi-square  value  which  is  a measure  of  the  degree  of  consistency  between  the 
measurement  model,  the  input  estimates,  and  the  covariances  of  the  input  quant- 
ities. 

The  adjusted  estimate  of  the  value  of  a measured  quantity  differs  from  the  input 
estimate  only  if  the  measurement  model  imposes  additional  information  regarding  the 
value  of  that  particular  quantity.  In  that  case  the  standard  uncertainty  of  the  adjusted 
estimate  will  be  smaller  than  the  standard  uncertainty  of  the  input  estimate.  For  a 
good  measurement,  the  difference  between  the  adjusted  estimate  and  the  input  estimate 
of  a measured  quantity  should  not  be  large  compared  to  the  standard  uncertainty  of 
that  difference.  It  has  therefore  been  suggested  that  the  ratio  d of  the  difference  to  its 
standard  uncertainty  is  calculated  and  assessed  against  a selected  criterion,  e.g.  \d\  < 2. 
By  plotting  the  d values  of  the  adjusted  estimates  it  is  possible  to  assess  whether  a too 
high  chi-square  value  is  caused  by  a few  poor  input  estimates  or  is  due  to  a poor  model. 
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